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I. INTRODUCTION 


This report describes the research performed under NASA Grant 
No. NSG 1004. The research utilizes an integro-differential method 
for numerically solving unsteady incompressible viscous flow problems. 
The method was developed previously by the principal investigator and 
his co-workers at the Georgia Institute of Technology. During the 
course of this research, several important advancements were made and 
incorporated into the solution procedure. Most of these refinements 
are described in arti::les that appeared in the open literature (Referenc 
1 to 4 ). A computer program was prepared using the integro-differential 
method to solve a problem of an impulsively started airfoil. The 
results for the airfoil problem have not yet appeared in the open 
literature. This report therefore consists of two parts. In the first 
part, the work that has been documented in the open literature are 
summarized. In the second part, some of the results obtained for the 
airfoil problem are discussed and ccmipared with available information. 
The second part of this report on the airfoil problem will be supple- 
mented by a Ph.D. thesis being prepared by the second author of this 
report and by journal articles. 
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II. THE INTEffilO-DIFFEREKTIAL METHDD 

The follotrlng articles were published a^ results of research 
supported by the present grant: 

(a) "A Plowfield Segmentation Method for the Numerical Solution 
of Viscous Flow Problems” by J. C. Wu, A. H. Spring emd N. L. Scudcar, 
Proceedings of the Fourth International Conference on Ihimerical ^thc^ 
in Fluid Dynamics, pp. U52-U57, Springer -Verlag, 1975. 

(b) "Integral Representation of Field Variables for the Finite 
Element Solution of Viscous Flow Problems," Proceedings of the 1974 
International Conference on Finite Element Methods in Engineering, 

pp. 827-840, Clarendon Press, 1974, by J. C. Wu. 

(c) "Velocity and Extraneous Boundary Conditions of Viscous 
Flow Problems," AIAA paper 75-47» 1975 j H pages, by J. C. Wu. 

In the remainder of this part of this report, the integro- 
differential formulation and the advancements reported in the above listed 
articles are summarized. 

The Integro -Differential Formulation 

The Navier-Stokes equations for a fluid with constant density 
p and viscosity v, and subject to negligible bcxiy forces, are expressible 
in terms of the velocity v and the pressure p as 

^ + (v . v)vs--VP + Wv 

Oc p 


( 1 ) 
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Together with the continuity equation 

? • V * 0 (2) 

the equations form a set of four scalar differential equations which are, 
in principle, sufficient for the determination of v and p, provided that 
adequate initial and bouMary conditions are known. Alternative systems 
of differential equations can be written in terms of the velocity and the 
vorticity, the stream function and the vorticity, or some other dependent 
variables, with each of these systems, it has been necessary to con^te 
the values of the dependent variables, usually by means of a finite 
difference method, for the entire flow field. Associated with this common 
feature for previous methods there are at least two difficulties in 
treating problems involving the flow past the exterior of a finite body 
such as an airfoil, particularly for cases where the shape of the body 
is not geometrically simple. 

The first difficulty arises because the flow region is infinite 
in extent (at least in the mathematical sense). Consequently, boundary 
conditions imposed at infinity need to be satisfied. The second, and 
more serious, difficulty arises bacause of the large number of discrete 
data points required in the numerical formulation of complex problems. 

In particular, for the general problem of flow past a finite body, there 
are a boundary layer and a sepeirated flow region both embedded in a 
potential flow. The magnitudes of the velocity gradients differ greatly 
in these different flow regions. To obtain a sufficiently accurate 
resolution of the numeric.il solution in all three flow regions, a 
large number of discrete grid points must be employed. The data storage 
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and counter time requireronts are, therefore, also leurge. 

Zt is convenient to introduce, as a dependent variable in the 
problem, the vorticity vector defined by 

0) « 7 X V f3) 

and consider the vorticity transport equation 

^ » 7 X (v X 3) + vV^3 (4) 

which is obtained from equation (l) and describes the convection and 
diffusion of vorticity with time, t. 

The set of equations (2), (3), and (4) partitions the problem 
conveniently into a kinetic part and a kinematic part. 

The kinetic part of the problem deals with the change of the 
vorticity field 3 with time and is described by equation (4). Since 
equation (4) is parabolic, with known spatial dlstributicxis of v and 3 
at any given time t , the rate of change of 3 with time and h-ence the 
new distribution of 3 at a subsequent time, say t + ^r» can be confuted. 

It can be easily seen from equation (4) tlmt, in the computation of the 
new vorticity values, it is necessary to know tte values of velocity at 
the given time t onl^ in the region of non-zero vorticity. Since this 
region of ron-zero vorticity is the only region where the flow is 
rotational and hence the viscous forces are non-aero, the kinetic part 
of the computation can be confined to the viscous region only. 

The kinematic part of the problem relates the velocity distribution 
v(r, t) at any given instant to the vorticity distribution 3(r> t) at 
that instant and the velocity cor^itions at the boundary of the flow 
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regicm, This kin^aatlc relation consists of equations (2) and (3) 
together with prescribed velocity boundary conditicm. It Is esqpressible 
as a vector R>isson*s equation for the velocity vector 

q^v(r, t) o - 7 X 3(r, t) (5) 

together with the prescribed boundary conditions for v. Equation (5) 
is obtained by taking the curl of equation (3) and using equation (2) 
to simplify the resulting equation. The vector r is tte position vector. 
Since equation (5) is elliptic, the con5»utation of v at the neighboring 
data points. It is, therefore, not possible to coiapute v ea<plicitly 
(point by point) in the flow field. That is, if equation (5) is xised 
to compute the values of v frcm known values of o) and known bouz»iary 
conditions, the computation must be performed for the entire flow field, 
inclusive of the viscous region where 3 is non-zero and the inviscid 
region where 3 is zero. Using an integrstl representation of the velo- 
city vector, however, it becomes possible to evalimte v explicitly, 
point by point. The need to compute v in the zero 3 region is eliminated. 

In the eeoplier developnent of the integro-differential method 
(Reference 6), it has been shown that, for the external flow problem, 
the integral representation of v is 


7(r, t) - - i J 


S(? . t) X (r„ - r) 


I? - ?1'^ 

I o ' 


dR + V 
O to 


( 6 ) 


where in three-dimensional problems A *= 4 tt and d = 3> in two-dimensional 
problens A 2n and d * 2, in one -dimensional problems A « 1 and d « 1; 

V is the free stream velocity; R is the enuire region occupied by the 
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fluid. The no->slip condition, v « 0 on the solid surface, and the velo- 
city boundary condition infinitely far away from the surface are incor- 
porated into the integral representation (6). 

With this integral representation, v(r, t) cam he computed 
e^qplicitly, i.e., the evaluation of v at each data point is accomplished 
independently of the evaluation of v at any other point in the flow 
field, by numerical quadrature of the integral in equation (6). As a 
consequence, in the kinematic part of the problou, the computation of 
V can be confined to the viscous region only. 

Since, with the integral representation, both the dynamic and 
the kinosatic parts of tl% problem can be confined to the viscous region 
of the flow, the nxsaber of data points required by the integro-differen- 
tial method Is drastically ssmller than that required by current o»thods 
treating the kinematic equation in the form of Poisson’s equation. For 
external flow problems, the vlscoiis region usually occupies only a 
small portion of the entire flow field and the limprovanent in compu- 
tational efficiency is particularly dramatic. For exanple, numerical 
solutions for a highly complex three-dimensional problem of a jet 
issuing from a plate and interacting with a crosswind over the plate 
were successfully obtained vising the integro-differential method 
(Reference 6). 

For the external flow problem associated with an immersed solid 
body whose velocity rises suddenly at an initial, time from zero to a 
constant value, the initial velocity field is given by the potential 
flow solution everywhere except on the solid surface. Tte no-slip 
condition at the solid surface gives rise to an initially concentrated 
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sheet of vorticity distributed on the surface, the vorticity everywhere 
else being zero. As time progresses, the vorticity diffuses outward 
fr<»ii the surface Into the fluid and subsequently also convects with 
fluid motion. The Integro -differential method allows the time develop- 
ment of both the vorticity and velocity distributions to be followed 
computationally, offering. In addition to the time -dependent solution, 
steady state or periodic (KdmuCn vortex shedding) solutions In the limit 
of large time (Reference 5)- 

The ability of the integro-differentisd method to confine tl^ 
st\;^y of a given problem to the viscous region is distinct frcas the 
basic notions of the boundary layer theory, or of the more general 
matched asymptotic e^^panslon method. With the Integro-dlfferentlal 
method, the viscous euid the Invlscid regions are not studied separately’. 
Neither a "matching" nor a "limiting" process is used. The ability to 
confine the study to the viscous region resul-ts not Ttobl a simplifi- 
cation of the Navler-Stol^s equations but l^om the use of the integral 
.representation. Since no simplification or approximations other than 
those Involved in the derivation of the Navier-Stokes equations in 
tlmlr differential form are introduced in the Integro-diffez^ntial 
formulation, the method is valid whether or not appreciable flow 
separation occurs. Furthermore, no information is lost by confining 
the computation -bo the viscous region. The problem of satisfying the 
boundary coniition imposed at infinity is removed. Information about 
the inviscid region is implicitly contained in the knowledge of the 
vorticity distribution, which is non-zero only in the viscous region. 
That is, the velocity in the inviscid region can be ccmputed if desired. 
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for equation (6) ia valid everywhere in th/6 flow field , altlKJugh in 
following the development of the flow field with tine such conputationa 
are unnecessary. 

The Flowfield Segmentation Method 

The earlier approach of the integro -differential method outlined 
in the proceeding section permits the computational field to be confined 
to the viscous region. The resulting reduction in the needed number of 
data points makes the integro -differential method much more efficient 
than other methods. More recently, it has been shown that a generalized 
integral representation of the velocity vector further permits the 
viscous region to be divided into segs^nts of arbitrary shape cuid size, 
and the conputation to be performed within each segment individually. 
Error cuialyses and ninmerical illustrations show that this flow field 
segmentation tecMlque leads to further substantied Improvoiient in 
solution speed and accuracy. 

The generalized integral representation is derived from equations 
(2) and (3) by using a vector potential. Green’s theorem for vectors, 
and principal solution of the Laplace's equation and is 


▼ (?) - - r I 1 — ^-3— ^ — T 

* L Jr . ?|'l o Jr |r. - rl* ° 


(7,xn^)x(?^.?) 

l?o - 'I* ° 


( 7 ) 
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vhere the sutnicript denotes that the variable is in the space. 

This representation is applicable to the internal flour problem 
as well as to the external flow problem. For the external flow problem^ 
if one c<msiders R to be the entire region occupied by the fluid) then 
the integrals over B give - v^ eur^ equation (7) reduces to equation (6). 
Alternatively) one may consider R to be any portion of tl^ fluid domain. 
Then, Knowing tiie vorticity distribution in that portion of the fluid 
domain and the velocity on its bo\uidary) ths velocity distribution 
within that portion can be computed by using equation (7). 

This new representation therefore permi'to the flow field seg- 
mentation for tt]» external flow problem as follows: 

In the kinetic part of the computation) knowing the distributioxu 
of V and m at points interior of and on the boundary of a segn»nt at a 
given tjjae, new vorticity values within the segment at a subsequent 
time are counted by numerical solution of the vorticity transport 
equation (U). Computation of xmw vorticity values at the boundeury 
points are obtained by either a numerical interpolation) or an explicit 
finite difference method) or by treating tl^se points as interior points 
of a different set of segments. Exceptions to the above are the boundary 
points on solid surfaces or openings where specified velocity boundary 
conditions enter the problem and only the interpolation method can be 
used. The degree of accuracy of the boundary-points computation is 
made compatible with that of the interior point ccxnputation. 

In the kinematic part of the ccxoputation, new values of the 
velocity at the boundaxy points are computed from the new vorticity 
distribution explicitly by numerical quadrature of the integral 
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representation (6). With the new velocity at boundary points and new 
vorticity at the interior and boundary points determined, the computation 
of the new velocity values for ths interior poixits can be perform«i 
either by the numericeLl quadrature of the integral representation (7) 
or by a numerical solution of the Foisson*s equation. 

The advantage of employing the segasentation method can be easily 
seen frem the different numbers of algebraic operations x»eded fc»r 
computing each yelocity^ value using the integral representation. For 
illustration, consider the two-dimens i<mal flow problem, fiy mapping 
the solution field into finite-elements, one obtains the following 
formula for nurorical quadrature from equation (6): 


lA 


N 

I 

J»1 


F 


t,id 


(8) 


where denotes the 4th component of v and is either 1 or 2, ”i” 

refers to the node for which the valv^ of v^ is to be c<spu^d; 

refers to any non- zero vcnrtlcity node; F. . . are geometric functions 
depending on the relative positions of the nodes i and N is tl% 
total number of non-zero vorticity nodes in the flowfield R. Similarly, 
one obtains from equation (7): 

Q 2 P 

“ Z ^lAS A Z ®m,ik\,k 

b,b 1 hpl 

where Gr .. are geoiMtric functions depending on the relative poslticms 
m,ih 

of the node i and the node k which U located on the compartment boundary; 
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Q is tbs number of mm-zero vorticity in t!«e region R| and P is the 
nuaber of boundary nodes on B. If tlusre are 4000 iK>deB on which the 
values of vorticity and Its derivatives are non-negligible, then velo- 
city cboiponents need to be computed only for tiMse 4000 nodes . Without 
the segmentation of the solution field, one would usr equation (8) to 
conqpute the velocity values. Since N ■ 4C00, the computation of each 
velocil^ value requires the multiplication of 4000 values of the 
geoQMtiric function F. . . with 4000 values of the vorticity and the 
addition of the products. For the 4000 i^es (and therefore 8000 
velocity veQ.ues), a total of 3*2 x 10^ multiplications and 3*2 x 10^ 
additions are needed for the velocity coaptation. With the solution 
field segmented, one first uses equation (8) to compute thr. velocity 
values at nodes on the boundary of the comparteients . This portion of 
the velocity computation still requires 4000 multiplications and 40CX) 
additions for eacn velocity value. However, for the coaptation of 
velocity values at nodes interior nodes, equation (9) ia now used aikt 
the required number of su^tiplications and additions reduces to Q 2P 
for each velocity value. For exBaQ>le, with a compartment containing a 
total of 400 nodes, 80 of which are on the boundary of the conprtments, 
only 560 multiplications and additions are now required for coopting 
each /elocity value at an Interior node. Obviously, the segmentati^ 
of the solution field offers a very drastic ia^rovement in solution 
sy-eed. This I mp rovement can, in fact, be made even greater by the use 
of successive segmentation. That is, one may first coepte the velocity 
values on boundaries of larger conprbaents each coi^alning, say, 2000 
nodes. Once this is accomplished, the couptation of the velocity 
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values on the boundaries of the 40O-node conSMurtronts (but are interior 
of the 2000-node coBqpartments ) can be acconplished much more rapidly. 

The flow field segmentation technique utilizes the most important 
feature of the integro-differential method, i.e., the ability to confine 
the computational field to the viscous flow region only, and adds a 
great flexibility to the method. Since the segments can be of arbitrary 
shi^ and size, single shi^s (e.g., rectangular, triangular, circular 
shapes) can be utilised. Consequently, nMthods which possess great 
speed and accuracy but are difficult to apply in problepas of complex 
solid ociy shapes can be readily incorporated in the coimputati«Ml 
procedures. For exaa^le, the fast Poisson solvers can be used In the 
proposed work In the kincmiatlc ccmputatlon of the velocity at tim 
interior point of each segment, liie alternating direction implicit 
method can be used for tbs dynamic ccmiputation of the vorticity* 
Furthexmore, since the computation is performed likLividually within 
each s^ment, the use of variable grid spacing is straightf<»rward. 

Thus, for the segswnts adjacent to tte solid body ^diere the velocity 
and vorticity gradients are expected to be large, closely spaced grid 
points can be used to obtain sufficient solution resolution. 

Ihanerical results obtaii^ utilizing the segn^ntatlon method 
are found to be practically identical to those obtained without using 
segmentation. The results, discussed In more detail in Ref. 1, also 
show that very drastic reduction in needed computer tii)» resulted from 
the use of the segmentation technique, as was expected. 
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FarfieldL Bomatoy Condition 

It isdiown in ReferexK^e 3 that, for the external flow problaa, 

the difficulty of satisfying the freestream velocity boundary condition 

at an infinite distance from the solid body is ccmpletely removed by 

the use of the Integro-differential formulation. If the solid body Is 

not rotating, but is merely translating at a velocity - v (which may 

€0 

be timeHlependent) relative to the freestream, then eqnntion (6) Is 
applicable. The contribution of the freestream velocity boundary 
eolation to the velocity field is simply v . The vorticity 3 in tt^ 

CO 

fluid originates from the solid surface and, at any finite time t, is 
non-negligible only within a finite distance from the solid surface . 

Thus the contribution of the integral in equation (6) to the velocity 
field is zero at infinity. That is, with equation (6), the freestream 
velocity boundary condition is satisfied exactly at an infinite distance 
from the solid body. 

In Reference 3> the importance of the ”-**arfield" velocity 
boundary condition is illustrated by numerical examples. In addition, 
integral expression for the velocity vector are derived for the more 
general case where the solid body is undergoing both translation and 
rotation. It is shown in Reference 3 that the "farfield" velocity 
boundary condition is satisfied exactly by the integral es^ressions 
for this more general case. 

Extraneous Boundary Condition 


The kinetic processes of vorticity diffusion and convection in 
an incompressible fluid is described by equation (4), the vorticity 
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transport equation, ^is equation is parabolic in its tliae-space 
relation and is valid in the fluid domain H. In order to solve equation 
(4) and obtain the time dependent vorticity field i^(r, t) in the fluid 
domain R, the values of m on the boundary B must be known as a function 
of time. The boundary B may contain a solid boundary S. The "boundary 
values" of m on S may not be computed by solving equation (4) since the 
locid "generation" or "depletion" of vorticity on S is not described 
by the kinetic processes of vorticity transport. Since the specification 
of vorticity values on S apparently overspecifies the problem, the 
e:q>ression "extraneous boundary condition" is used to describe this 
difficulty. 

In most of the previous investigations, various one-sided 
difference formulae have been used to estimate the vorticity values on 
solid boundaries from the no-slip condition and the computed velocity 
(or stream function) values at grid points near the boundary. It has 
been found that difference formulae of first-order accuracy tend to 
yield stable solutions while second-order formulae tend to give unstable 
results at high Reynolds numbers. In fact, for some problems, second 
order formulae, even when stable, gave less accurate solutions than 
did the first-order formulae. The use of first-order formulae, however, 
restricts the over-all accuracy of the solution. There existed, there- 
fore, considerable uncertainties regarding the correct one-sided 
formulae to use and the limitations of each formula. 

With the integral representation for the velocity vector, it 
becomes possible to ccanpute tte vorticity values on S accurately based 
on kinematic considerations alone. In Reference 3> it is demonstrated 



that this new approach correctly simulates the phi^icfiJ. process of 
vorticity "generation** on the solid surface S and therefore removes 
the difficulty of extraneous boundary condition. 

The kinosatic aspect of the viscous flow problem relates the 
velocity distribution at any given instant of time to tl» vorticity 
distribution at that instant, and vice versa. If the velocity field 
is known, then the vorticity field is uniquely determined on the basis 
of equation (3). On the other hand, if the vorticity field is known 
throughout an unbounded xegion, then the velocity field is uniquety 
determined on the basis of equation (6). It is obvious that the stress- 
strain relation, which differentiates the kinetic behavior of a solid 
from that of a fluid, does not enter into the kinematic relatioiwhip 
between the instantaneous vorticity and velocity fields. Therefore, 
equation (6) is valid in both the fluid dctmin and the domain cx:cupied 
by the solid. Ih other words, as far as the kinematics of the problem 
is concerned, one may replace the solid by a fluid laovlng at the 
Instantaneous velocity of the solid. In a reference frame attached to 
the solid, the velocity inside the region occupied by the solid hapx>ens 
to be zero. The vorticity in the solid region is therefore also zero. 

As a consequence, the velocity field is uniquely determined by equation 
(6) both in the fluid domain R and in the solid dcanain R* provided that 
the vorticity field is known throughout the fluid dcanain R. The kine- 
matic relation requires that this vorticity field must be such that 
the velocity field as determined by equation (6) is zero everywhere 
inside the solid domain R* . Clearly, not all vorticity fields in R 
would give zero velocity values throughout R’. In other words, the 
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vorticlty field in R is subject to a kinematic restriction. 

One notes that, in addition to this klneoatic restriction* the 
vorticity field in the fluid domain R, except at the fluid-solid inter- 
face is determined by the kinetic process of vorticity transport. 

As a consequence* the fact that the velocity field is zero in the solid 
domain R* pinces a kinematic restriction only on the >;ortlcity distri- 
bution at tlw surface S . This restriction* in fact* gives rise to the 
phencmiena coomonly referred to as the generation or depletion of vorti- 
city on the solid surface. 

In Reference 3» it is fouM that for two-dimensional flows* If 
a vorticity field in R is found such that either the normal ccmiponent 
or the tangential conqponent of the velocity is zero on the solid surface 
S* then both velocity ccaaponents are zero in the entire region R* bounded 
externally by S. As a consequence, to establish the extraneous vorticity 
boimdary condition on s’*", one needs only to search for a vorticity dis- 
tribution on which* when combined with the known (from the kinetics) 
vorticity field everywhere in R except S , gives either zero nom»l 
velocity on S or zero tangent iaJL velocity on S. In other words, for 
the external flow problem, equation (6) is rewritten as 


R-S 


5 X (r - r ) 


dR. 


- ~ r 

2tt J . 


Cq X - r) 
I"* 

l'■o - ■'1 


dS + V 

O 00 


( 10 ) 


where ^ is a surface distribution of vorticity on S , treated separately 
from the remaining vorticity field in the region R - S . The vector 
equation (lO) gives two scalar component equations. With the velocity 
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bouzKlary ccmdltlon specified S and the vorticity distribution in the 
region R - known, the application of equation (10 ) at each bouMary 
point gives two equations, one for the normal velocity and one for the 
tangential velocity, both containing the values of Q on the boundary 
points as unknowns. If N boundary points are used in the numerical 
procedure, then the application of equation (10 ) at the N points gives 
a set of 2H algebraic equations containing N unknown values of Q. 

However, since one needs only to satisfy either the nom^ velocity 
boundary condition or the tangential velocity boundary condition, only 
N algebraic equations are needed in the computation. 

It is further shown that for an internal flow problem, the 
solution of the N algebraic equations exists and is unique. For an 
external flow problem, however, the solution of the N algebraic equations 
is not unique. The coefficient matrix of the N equations is in fact of 
rank N > 1. It is shown that for this external flow probl«n a unique 
distribution of the vorticity on the solid boundary is obtainable from 
the integral representation for the velocity vector and the integral 
law for vorticity. The analyses for the three-dimensional and interior 
flow problems are similar to those given here for the two-dimensional 
exterior problem. In fact, they are slightly less involved since for 
the two-dimensional exterior probloos the flowfield is doubly connected 
and for the other problems the flowfield of concern is simply connected. 
It is worthwhile noting that the well accepted Mthod for initiating 
the time -dependent viscous flow ccmputations is to assume that the 
solid is set into motion impulsively, obtain a potential flow solution, 
and compute the initial vorticity distribution (which is non-zero only 
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at the solid boundary) fraa. the tangential velocity of the potential 
flow on the boundary. The initial vorticity distribution is thus deter- 
mined from the kinematical aspect of the problem. This method for 
Initiating the ccm^tation is in fact a specialization of the present 
method to the case where the vorticity is non-zero only at the boundary. 

The present method does away with the need of using one-sided 
differeiKse formula. It allows the generation (or depletion) of vorti- 
city at the solid boundary to be followed ccm^tationally using kine- 
matic relations only. The present method ensures the satisfaction of 
the principle of total vorticity conservation auad of the specified 
velocity boundary coxviitions at all data points on the solid boundary. 
With a one-sided difference formula, the no-slip condition is applied 
ii^ividually at each grid point on the solid boundary and a vorticity 
value is obtained for that point. Each vorticity value thus computed, 
however, gives rise to a finite velocity at all the other boundary nodes. 
The sum of the contributions of boundary vorticities does not in general 
vanish at a^jr of the boundary nodes. Thus, in reality, there is no 
assurance that either the velocity boundary condition is satisfied or 
the total vorticity is conserved. 

It should be noted that Lighthill (Rference 8) described an 
approach which appears to be similar to that described here. In 
reality, however, there are substantial differences between Lighthill *s 
approach and the present work. These differences shall be e:q)lained 
in a future article. Thus far, with Lighthill *s approach, it was only 
possible to treat problems involving simple two-dimensional gecmietries 
such as circles and flat plates. The present approach has been used in 
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the study of problems Involving more complex geometries, including the 
three-dimensional problem of a jet issuing noxmally from a plate into 
a corssflov (Reference 6). Ifumerical illustrations are provided in 
Reference 3* 


Integral Reiafesentation for Ccmipreasible Flows 

In Reference 2, an integral representation is presented for the 
kinematics of confess ible flow. This integral representation is 


7(?, t) = - i [ J 
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where 


b = 7 • V 


is the dilatation of the fluid. 

Reference 2 describes a procedure for coopting time-dependent 
viscous compressible flaws using the integral representation (U). 

Integral Representation for the Kinetics of Steady Incompressible Flow 
In addition to the integral representation for the kinematics of 
compressible flow, it was found that for steady flows, the kinetic 
aspect of inconpressible flow can also be recast into the form of an 
integral represention — for the vorticity vector. The result is, 
for the steady incon^ressible flow problem, a system of integral equati(xis 



i^ally suited for the finite-el^i^t approach. Each inte^al eqiuation 
la an Integral representation of a field variable. In the derivation 
of the integral equations, neither the existence of a variational 
principle nor a weighted residual (or Gad.erkin*s) approach is required. 

Hie application of the finite eluent method Is tterefore indepei^ent 
of a minimization of functionals. This new fozmulation for the steady 
flow problem retains the distinguishing feature of the integro-differentia] 
formulation for time -dependent flows — the ccaQnitation field can still 
be confined to the viscous region only. A method for coi^puting steady 
flows based on the new formulation is described in Reference 2. Scow 
numericed results are presented. 
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III. FLOI^ HIST AH BIFUISIVELY SMARTED AIRFOIL 

This part of this report descrlt^s the work that has been carried 
out In obtaining a numerical solution of t^e incon^resslble viscous flow 
past an IcqpulslYely started 9^ thick synnetric Joukowsky airfoil at an 
angle of attack of 1^^ and a Reynolds number of 1000. This problem was 
studied previously by Mehta (Reference 7) ifho mapped the fluid dcmialn 
into the Interior of a unit circle and used the vortlcity and stream 
function as the primary variables . 

To facilitate a coiqparlson between the results of this research 
and the results of Metha, the Intcgreil representation and the vorticity 
transport equation presented earlier in terms of the velocity and the 
vorticity are modified and reformulated in terms of the vorticity and 
the stream function In Mehta's transformed plane. Results are obtained 
by numerically solving the resulting integro -differential equations for 
S(»« early time levels and compared with that of Reference 7. 


Ctoordinate Transformtlons 

Ihe shape of the 9^ synmetric Joukowsky airfoil Is given in 
Figure 1 in the conq^lex plane z » x -i- ly. Qy using the conformal trans- 
formation 


1 ^ ^ 
— + V + 
H ' 


1 + YH 


( 12 ) 


id 

where k « re , C and y a-re real constants, the airfoil surface trans- 
forms into a circle and the exterior of the airfoil, i.e., the fluid 
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domain, la ma^d In k plane aa the Interior of the circle. 

Ihe coordinatea of the a-plane and n-plane are related by 

•In 0 / C“r , N /, 1 . \ 

y - — ( ^ - 1 } (l>^) 

vhe]^ 


2 2 

1 + Y r + 2y r cos 6 


(15) 


The scale feustor H la gl\ren by 
Where 

■ 1 + Y^ + 2y r cos 0 (17) 

For the present probloa, the values of C and y ^re taken to be 
0.92hl63^ and - 0. 0^214 respectively. The coordinates and the flow 
variables are to be nonnliaenslonallzed with respect to the radius of 
the circle and the freestream velocity. Therefore the airfoil in the 
z plane transfonos Into a unit circle in the x plane. The non-dimensional 
chord length, L, of the airfoil Is, according to equation (12), 

3.71281277 . Ohe airfoil Is 9-999&H> thick. 

A non-conformal transformaticni is utilized to obtain an efficient 
distribution of grid points for ccmpitatlon. The 0 coordinate in the 
X plane Is not modified. The r coordinate was "stretched” by the relation 
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p - [ tanh"^ (1.996590918 r - 1.032563339) + 2.8 ] (18) 

This transfomation msps the annular region betwe«i r » 0.02 and r ■ I 
in the K-plane into the interior of a unit circle in the p-6 plane. 

In Reference 7> both the kinonatlc and kinetic eqtuations are 
expressed in their differential form in the p-6 plane. Finite difference 
methods are then used to obtain numerical solutions. In the present 
work) the kinetic equation, l.e., the vorticity transport equation, is 
^qpressed in the differential form in the p-6 plane. The kinematic 
equation, i.e., the Poisson's equation, is recast into an integral 
respresentatlon for tlM stream fuz^tion in the h plane, as discussed 
in the following section. 
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The Poi88on*e equation in the k plane is 

^ i ^ + 4 . Bl^ (20) 

8r " r** ®e 

In Reference 7> equation (20) is re«*e3qpressed. in the p-9 plane. 
"Farfield** boundary conditions for the stresm function is then used at 
the boui^ary p « 0^ which corrdsponds to r «■ 0.02 in the K-plane and 
roughly corresponds to a circle of radius 13L in the z plane (the 
physical, plane). In the present study an integral representation for 
^ is derived and the freestream condition is satisfied infinitely far 
frm the airfoil (i.e., at r » O). 

In general, the Poisson's equation can be recast into an integral 
representation for the depeMent variable by the use of the Green's 
theorem and a principal soluticm of the Laplace's equation (Reference 
?). For two-dimensional problems, the prisuiipal solution corresponds 
to a logarithmic potential (Befer^ce 9)* integral representation 
for the Poisson's equ%tion in the z plane, ■> - m» is 


♦ (*) ■ J J «o ' “o 

R 




( 21 ) 


where 


f 




( 22 ) 


is the principal solution in the z plane. 

Using the prescribed boundary coi^itions. 
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- 0 


oa tba toUd •urface S 


(23) 



lAflnit«ly sway froa 8 (24) 



it can be ehown that the second integral in equation (21) reduces to 
- v^x. nserefore one hu 

♦(») - J J »o ^ "o * V - V + h 

R 


^ere is an arbitrary constant. 

With any confomal transfbznation z ■ g(K)« the transfoxmed 
Poisson's equation is 

V^il * - 0^ (26) 

vbere H is the scale fact^ and the subscript k denotes differentiations 
in the H plane. 

Equation (26) is still a Foi/(son'b equation. It differs froia 
the Poisson's equation in the z plane only by the fact that the 
inhomogeneous texm in the transforsttd equation is - in place of 
- m. Thus, one may let 

( 27 ) 


and %rrite 
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-11% 'k «0 

R ® 


irticre 



( 28 ) 


f 

K 



1 

K " 


Mid tlM su^cript ”k” refers to the k plane, e.g., is the region in 
the K plane that corresponds to the region R in the z-plane. If the 
transfoxmticm is such that the scale factor approaches unity for large 
valAies of z, thMi tlM streaaliz»s at great distance fron S are not 
distorted by the transfon&ation. The surface integrals in equation 
(28) then give Ujn ~ ^re % and t) are respectively the ahscissa 
ai^ ordinate in the transfon&ed plM». In the xnresent study, the 
exterior of the airfoil is transformed into tte interior of a unit 
circle in the K-plaxie, Ihe contribution of the suz'face integral in 
equation (2$) im.\st tterefore be analyzed sei«rate3y. 

Consider the surface 8* to be a circle of radius A centered 
about the origin, with A Tbt surface S' encloses the airfoil in 

the r plane. The stream function on 8* is - Vjc. In the trans- 
fon^ plane k, a circle of radius s centered about tte origin, 

with c i/R 0. Since f is a scalar and is invariant to the trans- 
f<»rmation, one obtains txm equation (24), 



Together with the boundary conditions on S, the surf sice integral in 
equation (28) gives - l/U(U sin 0 + V cos e). Equation (28) thus 

• 00 CO 

becones 


R ° I* - 

K 

- - (U Sin e + V cos 0) + C, (30) 

r * o» X 

K 

The abritrary constant is the value of the stream function outside 
the unit circle in the h plane. It is also the value of the stream 
function inside the airfoil in the z plane. The value of can be 
set to be zero and equation (30) can be more conveniently re-written as 

R ^ lx - *01 

K 

- ^ sin 0 + cos 0) (31) 

for reasons to be e3q>lained shortly. 

Although it is possible to re-e3Q>ress equation (31) in the p-0 
plane, for computational purposes it is more convenient to use equation 
(31). In the present study, therefore, equation (3l) is used. 
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Extraneous Boimdary Condition 

The vortlcity distribution uj in the h plane is restricted kine- 
satically by the requirement that f « 0 exterior of the unit circle. 

(Note that this region corresponds to the domain of the solid in the 
z plane.) At any given instant of time, the vorticity distribution 
away from the fluid -solid interface is determined by the kinetics of 
the probloa. The kinematic restriction therefore applies only to the 
vorticity distribution on the fluid -solid interface at that instant of 
time. In Reference an Implicit method utilizing the integral repre- 

sentation (6) is presented. The method is specifically developed for 
the problem formulated in terms of the vorticity and the velocity. It 
requires the numerical solution of a system of N simultaneous algebraic 
equations involving unknown vorticity values at boundary points on tte 
interface. In the following, an eaqplicit method for the computation 
of the surface vorticity values is described. This new method is 
useful both for the vorticity-velocity formulation and for the vorticity- 
stream function formulation. 

The differential equations describing the kinematic relation 
between the vorticity and the stream function (or velocity) are linear. 
Therefore the method of superposition is useful in computing the 
extraneous boundary condition. If one separates the vorticity on the 
surface S fraii that in the region R but not on S, then equation ( 31 ) 

H H 

becomes 


^ I J V 


R -S 


K K 



+ 


-1 r 

2tt J 


S 

H 
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- — (U sin 6 + V cos e) 

r « ea 


(3a) 


where C is a surface distribution of vorticity, its dimension beii^ 
that of vorticity times length. Let us consider in two parts 

C, = + Ck2 (33) 

where C i is the contribution of the freestream condition to the surface 
vorticity cmd is the contribution of the vorticity away from to 
the surface vorticity. 

The surface vorticity distribution corresponding to a potential 
flow past a circular cylinder is readily available in standard texts. 
Ui>on conformally mapping the exterior of the circle into the interior 
of a unit circle, one finds, 


= 2(U„ sin 0 + cos 0) 


(34) 


The stream function associated with the surface vorticity 
distribution s^ given above, shall now be determined for the 
region r i 1. (Note that h = re , and that the region r i 1 corres- 

ponds to the domain occupied by the airfoil . ) From equation (32 ) , this 
stream function is 



sin 0 + V cos 6 ) In 
o « o 


5 ^ 1 

[ r*^ + 1 - 2r cos(0^ - 0)]* 



2 i 

+ a) In [r + 1 - 2r cos(0^ - 0)]^ d0^ 


( 35 ) 
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2 2 ? 

where Q « (U + V )* is the freestream velocity magnitude and a = 

CD to 

tan”^ V Integrating by parts | one then obtains 




TT C08(e^ + ff) sin(6^ - e) 

+ 1 - 2r cos(0^ - 0) 


d0. 


Letting 0 = 0^ - 0, one has 




IT-0 [cos 0 cos(0 + a) - sin 0 sin (0 + a)] sin 0 

2 ■■■ 

r + 1 - 2 r cos 0 


d0 


The odd part of the integrand above does not contribute to the 
integral. The even part of the integreind is periodic in 0 with the 
period 2 tt. One tterefore 1ms 


T - Or sin(0 + a) 

•^1 ” TT 


! 


2 sin^ 0 


o r + 1 - 2r cos 0 


d0 


^ <tr Sln . (e -f , g) r f" di r" ^ cos_2| jp 1 

TT L , 1 o— — a ^ j. T _ o*. o -• 


or + 1 - 2r cos 0 ^^o r + 1 - 2r cos 0 


EQ)licit expressions for the definite integrals in equation (36) are 
available in standard mathematical tables. For r > 1, one has 


I » — sin(0 + a) “ (U sin 0 + V cos 0) 
lr' r « » 


(37) 


Thus the contribution of C to the stream function value is the negative 

k1 

of the freestream contribution to the stream function (the last part of 


p2TT 

equation 32). It is easy to show that J = 

o 

not contribute to the total vorticity. 


0 so that C 1 does 
'*k1 
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Let dC^ be the contribution of the elementary vorticity 

tt) dR located at h to the surface vorticity distribution C ^y‘ 
«o HO o ^ 

the method of images, one finds 


^=*2 = 


(r - 1) (B dR 

' O ' KO KO 


2n[l + r ^ - 2r cos(e - 0)] 


(38) 


By the use of the method of residues, it can be shown that the 
stream function associated with d^^ is ex€u:tly opposite to that 
asscMsiated with m dR . Therefore the combined contributions of 

«P MO 

dC rt and of m dR to the stream function is zero inside the airfoil. 

H H 

Furthermore, 




(r - 1) (JD dR 


2tt[ 1 + r ^ - 2r cos(e^ - 6)] 
'■ o o o 


d0^ = - to dR (39) 

o XO KO ' 


dC o together with tu dR do not contribute to the total vorticity. 
•k2 ko ko 

Obviously, if one lets 


p 

I* inA (r-l)todR 

r [at - -l r r o ^ "ko ho 

■ J - a. 


K X 

then the combined contributions of C „ and u) distributed in R - S 

k k X 

to the stream function is zero in the region r > 1 (inside the airfoil). 
Furthermore, together with the vorticity away from S gives zero 
total vorticity. One therefore has 


n 

R -S 


(r - 1) uj dR 

O ' "XDKO 


[1 + r ^ - 2r cos(0 - 0)] 
o o o 


+ 2(U^ sin 6 + cos 0) (4o) 
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Clearly, the surface yorticlty dlstrlbuticni given by equation (^) 
satisfied both the Idnematic requirement that the stream function 
inside the airfoil is zero euid the principle of total vorticity con- 
servation. It pexmits the e:^liclt, point by point, confutation of 
the surface vorticity values. 


Farfield Velocity Boundary Condition 

It can be sbowa that for r < 1, equation (36) gives 

I, a r(U sin 6 + V cos 0) 

J. CB 00 

so that I^ is zero at the origin. Furthermore, the ccaabined contri- 
iHitions to the stream function of «tnd the vorticity distributed in 

R - S is zero at the origin. Since the origin in the h plane corres- 
ponds to a surface infinitely far from S in the physical plane z, the 
freestream condition ^ - l/r(U^ sin 0 + cos 0) (or f - Vjc 

in the z plane) is satisfied exactly at infinity in the z plane. 

Aerodynamic Loads 

On the surface S the velocity vector vanishes and the Navier- 
Stokes equations reduce to 

^ss-^7^ v« - (4l) 

The tangential component of this equation is, in the h plane. 


dp h dw L , Su) 
d0 “ R dr " R P dp 


for r 5= 1 


(42) 


The pressure is computed from 
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9 

p(e) = p(0) + I P' J ( ^ ) for r - 1 




The pressure coefficient is defined by 


c (e) = 

® ipCu + V 2) 


(44) 


The skin friction coefficient is computed from 


Cj.(e) - - 1 01 


(45) 


The force coefficients normal and tangential to the airfoil chords axe 
denoted by Cjj and C^. The lift, drag, and moment coefficients are denoted 
respectively by Cj^, Cp, and Cj^, and are computed frcm and values. 

numerical Brocedure 

The grid system in the p-0 plane is formed by using dp = lA7 
and d9 = "n/hO. These values are identical to those of Reference 7» 

As shown in Figure 2, the coordinates of the grid point i,j are 
p s= 1.0 - (j - l) dp and e = (i - 1) d9. That is, the unit circle, 
p s= 1, which corresponds to the airfoil surface in the physical, plane, 
is assigned the value ;J = 1, the point p = 0, which corresponds roughly 
to a circle of radius 13 L in the physical plane, is assigned the 
value j = 48, the leading edge line is assigned the vfi0.ue i = 4l, the 
trailing e<3ge line is assigned the values i = 1 or 8l, as shown in 
Figure 1. Values of u)> are ccmjputed at the grid 



points i, ;J. Values of however, are caaputed midway between the 
grid points i,j and i *»*!, J, as shown in Figure 2 and designated .. 

For the kinetic part of the ccanputation, the vorticity transport 
equation, (l9)» is solved using finite difference methods. Both the 
ADI (alternating direction implicit) method and the SLB (successive 
line relaxation) method have been studied. 

With the ADI method, central differences are used for the 
diffusion term. For the convection term, four different differencing 
schemes have been tried. These are: (l) central differences, (2) 

second upwind differences (Reference 10), (3) combination of the first 
two (Reference 11 ), and ( 4 ) selective central or upwind differences 
depending on the cell Reynolds number (Reference 12). Scheme (l) gives 
spurioxis values of vorticity near the outer edge of rotational region 
around 0 » Schemes (2) and (3) give unsatisfactory Cp distribution 
on the lower surface near the forward stagnation point. Scheme ( 4 ) 
is found to be better than the other three from these respects. 

Using Scheme ( 4 ), the finite difference form of the vorticity 
transport equation is 


n+4 n 
0)^ I - o». 
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( 46 ) 


and 


n+1 n-*i 
“ « 




L / , \2 1 /n+1 p n+1 . n+1 v 

77:^ K.j-i - 


(AP)‘ 


^ L 1 / n+i p n+i n+4 % 


• (\) 3 * ♦ ",> 5 ‘ 

u 

Here, C^ and C are the finite difference form of the convection terms. 

These terms are written out as explained below. 

Two cell Reynolds's numbers based on the cell lengths in the 

e and p direction are defined respectively as 


and 


(■*Cg>i3 = Y'l >siti,d E “ 

P ^ 


( 48 ) 


( 49 ) 
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If the vulue of R_ <2, central difference is need for the term C.; 

otherwise ujwiM difference is used. Similarly, central diffei^nce 

for the term C is used if R. <2; otherwise u|iwind difference is 

p % 

used. These differences are given below: 


\ <2: (C^)^ - (rp‘)i^ ^ ^ ^50) 


R a 2 : 


Rc <2: 
P 


«C 

P 






(51) 






" ^ 3^*3 2hp ^^♦e‘“^i,j-i " ^♦e^^^i.d+i ^ 






1 • 


(53) 


* ’^a‘’’a\/“i.a-i ■ "la^ A? ’ “ *9 ^ ° 


In the above set of relations, the superscripts for oj indicating the 
time step ere as determined by equation (46) ax^ (47)- As of now, the 
stream function and its derivatives are allowed to lag by one time 
step and so the s\:q>erscripts for and are always n, irrespective 
of what is ingplied by equations (46) and (''■7). 

The solution procedure is as follows . For each value of J , 
equation (^) gives a set of 80 algebraic equations Involving 80 uidoKnm 
values of vorticity at the 80 grid points ”i" on the grid lines and at 
the time level (n + i^). The coefficient nmtrices for these sets of 
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ftl^ebralc et^uations it cyclic tridiagonal. Recursion relations (Reference 
13) are utilized to solve the sets of algebraic equation for each 
pertinent j value, beginning vith J « 2 and ending with the first grid 
line where tl» vorticity is zero at each grid point ”1” at the tijne 
level n. 

Equation ( 4 ?) is then solved to obtain the vorticity vulues at 
the time level n *1' 1. For each vedue of i, equation ( 4 ?) gives a set 
of algebraic equations involving unknown vorticity values at ths grid 
points and at the time level (n l). Ihere are 80 sets of algebraic 
equations, cn?e set for each for each value of i. The coefficient 
matrices are regular tridiagonal and Thomas algorithm (Reference 10 ) 
is used. It should be noted that the solution requires the surface 
vorticity at the (n + l)th time level be known. The surface vorticity 
values at the (n + l)th time level, however, depends on the vorticity 
values away from the surface kinematlcaJLly, as discussed earlier. An 
iterative procedure is tlmrefore adopted with which the surface vorti- 
city valiies coogputed frcm the kinematic relation, equation (39) » based 
on non-surface vorticity values obtained for each iteration are used 
in the subsequent iteration. For the initial iteration, the surface 
vorticity values are estimated from the surface vorticity values at 
time levels (n-l) and n through an extrapolation procedure. 

For the successive line relaxation schmae, using central 
differences for the space derivatives and forward differences for the 
time derivative, the residue is written eis 


(Aes) 


n+l,k+l 

id 


n+l,k+l n+1, k 
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m. 


n 






[°X 


d+1 


n+l,k+l + o+l.*»l ) 


, o^l.l^l n+l»fc*'l w c - 2a?T^ - «f*’f . ) 

“aK.j-i ■ "l.J+1 ’ 3 *« ^-i>l 


n+1 n+1 


■♦• C 


u(*o 

^ *’1+1.3 


n+1 

*1+1.3 


- t' 


n 


*’1-1.3 


n+1 

*1-1.3 


• '’"I,,-, cir - <sr > ] - « 


n+l|k 


(54) 


Where k 1® the iteration counter. 


C 


1 


I. 

R 


(4P)^ ' 


B Zip ' ®3 



=4 




The iteration counter® for the vorticity value® on grid line® 
i + landi-lin equation (48) are epecified in accordance with the 
direction of sweep. The sweep® always initiate wi.h the grid line 
1 B 4l (the leading edge line) and terminate with the grid line i « 8l 
(or 1). Thus, for the upper surface of the airfoil, the sweep is in 
the it direction, i.e., i:4l(l) 81. The iteration counter for vorticity 
values is k + 1 on the grid line i - 1 and is k on the grid line i + 1. 

For the lower surface, the sweep is in the i-^ direction, i.e., ij4l(-l)l. 
The Iteration counter for vorticity values is k + 1 on the grid line 
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i 1 and is k on the grid line 1. Ihe order in vhich the tq^per and 
loirer surface sireeps are done is alternated for successive iterations. 
For odd numbered iterations, the sweep along the upper surface is 
followed by the sweep along the lower surface. For even nuo^red 
iterations, the sweep along the lower surface is undertaken before the 
sweep along the upper surface. 

With the above prcwedure, a system of algebraic equations with 
a tridiagonal coefficient matrix is obtained from equation ($4) for 
etch value of i. Ttese equations are solved for values of 
not on the airfoil surface. The corresponding values of surface vorti- 
city are obtained by solving equation ( 39 ) ^ values of m at aU 

grid points ure tl^n re-set using the equation 


n+l,k+l 

*i,s 


n+l,k 


+ 0(A«) 


n+l,kfl 

i.d 


( 55 ) 


where p is the relaxation parameter and is assigned the value of O.09. 

The Calculatims with successive line relaxation, were rei^ated 
using Arakawa difference for tlw convection terms and thi^e-point 
teckmurd difference for the time derivative. The results presented 
here are obtained using Arakawa and three-point bac3o#ard diffez^nces. 

The method of coaptation of loads is identical to -^hat of 
Reference 7. 


Discussion of Results 

Solutions using ADI nmthod have been obtaii^d up to a time level 
of t B 0.212 and using SLR method up to t “ 0.02. Figures 3 to 6 show 
the vorticity profiles at several chordwise stations obtained using 
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tha ADI Dethod and the SIB method. A con^rieon of these solutione 
with those of Reference 7 shows that the ADI solution agrees well with 
tte SIB solution at^ the solution given in inference 7 at staticms 
away from the trailing edge (Figures 3> ? and 8). At the trailing 

edge station (Figure 6 and 7) the ADI solution differs substantially 
from that of Reference 7* In particular* the sl<^ of tlM vorticity 
profile changes drastically and repeatedly near the surface (Figure 6) 
so that no meaningful computation of the normal vo.*ticity gradient at 
the trailing edge is possible. Ihis eratic behavior* in fact* extends 
to chordwise stations near the trailing edge. In addition to the 
erratic instantaneous behavior near the trailing edge described above* 
it is found (although not shown here) that the vorticity values near 
the trailing edge as obtained by the ADI method oscillate violently 
as time progresses. As a c(»sequence* no reliable surface pressure 
distributim* which must be counted fron the normal velocity gradients* 
can be established from the ADI solution near the trailing edge. 

Ihe vorticity profiles at the trailing edge as obtained from 
the SIB metlK>d are in good qualitative agreement with that of Reference 
7. Quantitatively* there aie noticeable differences between the two 
solutions n^ir the surface which ii^reases gradually with the pro- 
gression of time. These differences are not une3g|)ected since tlm 
method of determining the surface vorticity in the present study is 
substantially different from that of Reference 7. 

The surface pressure distribution for ^ b 0.003* O.OOS and 0.02 
are plotted in Fig\ures 9> H and 13 . For the ADI solution the nman of 
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clockwise and antl-clockwise Integration is plotted. In the case of SLR 

solution the results presented correspond to integration from the leading 

edge along the upper and lover surfaces. In addition, the values are 

plotted so as to make the leading edge same in all cases. The values 

of the constant added to achieve this, are indicated in the figures. 

Figures 10 and 12 give the C distribution close to the trailing edge. 

Values of results obtained by Mehta are not included in these figures, 

as they are not available in such detail in Reference 7. As could be 

seen frcmi Figures 9"13 j the ADI method gives rise to meaningless values 

of C near the trailing edge. The ADI solution gives C curves that are 

not closed. That is, after the ccwipletion of a circuit around the 

airfoil the starting value of C is not returned. This is due to the 

P 

erroneous values of the vorticity gradient at the surface as predicted 

by ADI method. The ADI solution deteriorates further with time and 

gives rise to negative lift as could be seen from Figure 13. 

The SLR results for are obtained by integrating from the 

leading edge along the upper and lower surface. This gives rise to 

two values of C for the trailing edge. We define a tP obtained by 
P P 

dividing the difference between these two values by the iMocimum 
absolute value of on the surface. A comparison of these values 
with those of Reference 7 is as indicated below. 


Time Block 

No. of Time Steps 

Present Study, SLOR 

Average 

Mehta 

0.0 -» 0.006 

6 

2.161 

1.468 

0.006 - 0.012 

3 

0.640 

1.593 

0.012 0.02 

2 

0.400 

1.082 



The variation of Cj^ with time, for SLR scheme, is given in Figure l4. 
Compared to results of Reference 7> these values show a smoother 
variation with time. 

Fr<mt the above discussion it is concluded that successive over- 
relaxation is to be prefered over ADI method. As of now, the exact 
cause of the failure of the ADI method is not known. Ikifortunately, 
the SOR method is more time consuming than ADI scheme. It may be 
possible to reduce the computer time for the SOR method by clwosing 
a set of optimum relaxation parameters which has a spacewise variation. 
These conclusions are based on a few numerical ei^eriments carried out 
in advancing the solution from t = 0.002 to t = 0.003 . The following 
table summarizes the findings. 



Scheme 


Ifo. of Iterations 

Confutation 
Time (sec ) 

1. 

ADI 


40 

81 

2. 

Point SOR 


125 

178 

3. 

Line SOR with p = 0.09 


74 

i49 

4. 

Line SOR with p^ « 1.5 

60 

127 


In the above experiments, the convergence tolerance on tO) has been 
kept same at 0.5. In the ADI scheme the iteration indicated is for 
the surface vorticity during r-direction sweep. It is apparent that 
successive line over -relaxation has a slight advantage over successive 
point over -relaxation. The relation for given in item 4 above is 
an ad hoc selection and has no theoretical basis. 

In continuing this work, a proper investigation of the relaxation 



factor for SLR method Is to be done. There Is also a possibility 
of combined usage of ADI and Sl^ schemes, with the successive over- 
relaxation scheme limited to a small region near the trailing edge. 
This possibility is still to be explored. 
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